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Abstract 


Approximation methods are often used in porous electrode models to eliminate the need to solve the local solid phase diffusion equation. These 
methods include Duhamel’s superposition method, a diffusion length method and a polynomial approximation method which have long been used 
in the literature. The pseudo steady state (PSS) method is a method that has been used recently to develop a solution to the diffusion equation in 
a spherical particle with time dependent boundary conditions, but the PSS method has not been used in a porous electrode model. These methods 
are compared to each other in a dimensionless analysis study, and they are used in a porous electrode model to predict the discharge curves for a 
LiCoQ, electrode. Simulation results presented here indicate that the PSS method or the high order polynomial method should be used in a porous 


electrode model to obtain accuracy and save computation time. 
© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


Porous electrode models have been used widely to study 
the cell performance of various lithium based battery systems 
[1-8]. These models include mass transport in the liquid phase, 
electronic conduction in the solid phase and liquid phase in 
macro-scale and lithium diffusion inside solid phase particles 
at a local position in the porous electrode. The lithium transport 
in the solid phase is usually described by Fick’s second law for 
a spherical particle: 
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= Ds 1 
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The porous electrode model solution requires solving this 
particle diffusion equation in the an extra pseudo-dimension r. 
Consequently, the number of unknown variables in the model 
is much larger, which requires more computation time. There- 
fore, approximation methods are often used in porous electrode 
models to eliminate the need for the time consuming calcula- 
tion of solid phase diffusion in the extra r dimension. These 
methods include a Duhamel’s superposition method [1,2], a dif- 
fusion length method [9,10] and a polynomial approximation 
method [11,12]. In addition, Liu [13] developed recently an 
analytical solution to the diffusion problem with time dependent 
boundary conditions based on Olcer’s pseudo steady state (PSS) 
approach [14,15]. To our knowledge, Liu’s method has not yet 
been used a porous electrode model. In this work, these approx- 
imation methods are compared to each other in a dimensionless 
analysis study, and they are used in a porous electrode model 
to simulate the discharge profiles of a LiCoOz electrode. The 
advantages and disadvantages of these methods are discussed 
based on the results from both studies. This work provides useful 
information to help one choose an approximation method for his 
specific application to maintain accuracy and save computation 
time. 
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Nomenclature 

Brug Bruggeman coefficient 

Cmax maximum Li* concentration in the particles 
(mol cm~?) 

Cs Li* concentration in the solid particles 
(mol cm~?) 

Cs volume averaged Li* concentration inside 
spherical particles (mol em7*) 

Cge Lit concentration at the surface of the solid 
particles (mol cm~?) 

C dimensionless Li* concentration 

De diffusion coefficient of the electrolyte (cm? s7!) 

D; solid phase Li* diffusion coefficient in the 
particles (cm? s7!) 

f mean molar salt activity coefficient 

F Faraday’s constant (96,487 C mol`!) 

Jn pore wall flux (mol (cm? s)~!) 

k summation counter in Eq. (2b) represents 
previous time steps 

kp kinetic rate constant 
((mol (cm? s)~!) ((mol cm73)~ 1-5) 

ls diffusion length of a spherical particle (cm) 

n summation counter in Eq. (2b) represents current 
time step 

dm summation term used in (Eq. (9c)) 

a volume averaged concentration flux (mol cm~*) 

R gas constant (8.3145 J (mol K)~!) 

Rs radius of the spherical particles (cm) 

Rf contact resistance (Q cm?) 

S geometric area of the electrode (cm?) 

t time (s) 

i transference number of the electrolyte 

T temperature (K) 

Ueq equilibrium potential of the electrode (V) 

v thermodynamic factor of the electrolyte 

W active material loading in the electrode (g) 

x dimensionless distance in spherical particles 


Greek letters 


Qa, Ae 
B 
ô 
ôp» bs 
Ee, Es 


de, ds 
Ke 

Am 

Os 

T 


transfer coefficients 

symmetry factor 

dimensionless pore wall flux 

electrode or separator thickness (cm) 

volume fraction of the electrolyte or active mate- 
rial in solid phase 

liquid or solid phase potential (V) 
conductivity of the electrolyte (S cm7!) 
positive eigen-value determined from Eq. (9d) 
conductivity of the solid phase (S cm7!) 
dimensionless time 


2. Overview on approximation methods 
2.1. Duhamel’s superposition method 


Duhamel’s superposition method was the first approximation 
method [1,2] used in the porous electrode model. Duhamel’s 
method relates the solution of a boundary value problem with 
time dependent boundary conditions to the solution of a similar 
problem with time-independent boundary conditions by means 
of a simple relation [16]. The Duhamel superposition equation 
for the solid phase diffusion equation is [1,2]: 


9 
Dide) _ NECIT CoH) gy Com = Cond 4, 
Rs or r=R; k=0 At At 
(2a) 

where 
An-k = al(n — k)At] — a[(n — k — 1)At] (2b) 
and 
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An expression for a(t) was developed by using the Laplace trans- 
formation technique for short times and long times [1,2] and can 
be written in terms of a dimensionless time t = tDs/R?: 
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The expressions (Eqs. (3a) and (3b)) are both uniformly valid. 
However, Eq. (3a) converges more quickly with fewer terms at 
short times than Eq. (3b). 

The relation between the surface concentration at current time 
step Cs,n and time dependent pore wall flux jn can be established 
through Eqs. (1c) and (2a) as follow: 

Jn = ow m csk) 4 + (Cs,n = 


Cs n— 
R At nak o Pai 
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(4a) 
k=0 

It should be noted that values of the surface concentration from 
all previous time steps cs, are required in Eq. (4a) to calculate 
the value at current time step cs». The volume averaged concen- 
tration ¢, in the spherical particles can be readily calculated as 
follow: 


t 

3 

aac | -2 jy dt (4b) 
0 Rs 


2.2. Diffusion length method 


Wang et al. [9] used diffusion length concept to simplify the 
solid phase diffusion equation in the porous electrode model. By 
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assuming a parabolic concentration profile in the diffusion layer 
and using the volume average technique, they determined the 
diffusion length to be R,/5 for spherical particles. The solution 
of Eqs. (1a)—(1c) can be approximated with a set of differential 
and algebraic equations: 


OCs 3 


Be ps ea 5 
at R, Jn (5a) 

Ds — 
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where l; is the diffusion length and takes the value of R,/5 for 
spherical particles. 

The diffusion length method predicts that the surface con- 
centration and volume averaged concentration inside a particle 
are linearly dependent on each other, which should be valid 
only after the diffusion layer builds up to its steady state. There- 
fore, the method is inadequate at short times or under dynamic 
operations, such as pulse or current interrupt operations. 

In view of the shortcoming of the diffusion length method, 
Wang and Srinivasan [10] corrected the diffusion length method 
by empirically incorporating an intuitively expressed time 
dependent term into the diffusion length equations: 


OCs 3. 

“x R” (6a) 
D; 4 4/ Dst 

$ (cse — G3) = — jn { 1 — exp à (6b) 
ls 3 k 


The empirical term was formed based on the observation that the 
surface concentration increases exponentially at short times. The 
value of the multiplier in the exponential term was determined by 
matching the surface concentration profiles obtained with Eqs. 
(6a) and (6b) at given pore wall fluxes to those from Duhamel’s 
method. The value of 4/3 was found to be able to provide good 
results under a wide range of operating conditions [10]. 


2.3. Polynomial approximation method 


Polynomial approximation method [11,12] by Subramanian 
et al. (see also Rice and Do [17], Chapter 12) was also based 
on the parabolic concentration profile assumption and volume 
averaging technique. The authors first assumed that the concen- 
tration profile could be described by a second order polynomial 
(c=a+ br’) and the equations for low order polynomial method 
were derived as: 


OCs 3 
eee oe ey p! 
at R” (7a) 
5D . 

à (Cse — Ts) = Jn (7b) 
R; 


To achieve better accuracy at short times, the authors [11,12] 
used a higher order polynomial (c =a + br? + dr*) model for the 
concentration profile. The equations for this high order polyno- 
mial method are as follows: 

OCs 3 
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where q; is the volume averaged concentration flux which phys- 
ically defines the change of concentration with respect to the 
position in the system [12]. 

This high order polynomial method uses a different approach 
from the diffusion length method to improve the solution accu- 
racy at short times. The diffusion length method uses the 
empirical exponential term in the equation and determines the 
multiplier value by matching surface concentration profiles to 
the exact solutions. The high order polynomial method uses 
a higher order polynomial for the concentration profile in the 
derivation, and one could derive new sets of equations with an 
even higher order polynomial model, if needed, following the 
same procedures discussed in the papers [11,12]. 


2.4. Pseudo steady state method 


Liu [13] applied Olcer’s pseudo steady state (PSS) method 
[14,15], which is a form of a finite integral transform technique 
to eliminate the independent spatial variable r from the solid 
phase diffusion equation. For diffusion problems with a time 
dependent pore wall flux jn as a boundary condition described 
by Eqs. (1a)-(1c), the generalized PSS solution was found to be 
[13]: 


OCs 3 
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at R, Jn (9a) 
dqm Mn Ds 42 Dgt/R2 ; 
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z 2 
x Lin — eA Dst/R; qm] (9c) 


where qm is denoted only as a summation term in Eq. (9c) and 
itself does not have a physical meaning. The Àm values were 
determined from the eigenvalue equation: 


Am = tan(Am) (9d) 


Eqs. (9a)—(9d) were developed and implemented only for diffu- 
sion in spherical particles. However, the PSS method could be 
extended to diffusion in particles with other geometries, such 
as disks and cylinders, following the procedures detailed in the 
Ref. [13]. 


3. Dimensionless analysis 


Before these approximation methods are used in a porous 
electrode model, we use dimensionless analysis to compare 
how these methods perform on diffusion problems with time 
dependent boundary conditions. Eqs. (1a)—(1c) are converted to 
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dimensionless form using following dimensionless variables: 


Dgt 
Ca", #25, @ 2, sec 
Cmax Rs Re Dscmax 
(10) 


The governing equation and boundary conditions in dimension- 
less forms are: 


aC 1 2 3C 

E E saa 11 
ər x (x =) (la) 
ae) o (11b) 
əx | —0 

OC) Sio (11c) 
ax x=1 


The exact solutions to the dimensionless diffusion equations 
are compared to those obtained from approximation methods 
in Figs. | and 2 where the dimensionless surface concentra- 
tion minus the volume averaged concentration is plotted against 
dimensionless time. The dimensionless pore wall flux ô(t) takes 
the form of 1 — t and sin(5r), respectively, in the figures. The 
finite difference method with 100 internal node points was used 
to find a numerical solution to Eqs. (1 1a)—(1 1c) and is named the 
exact solution. Two summation terms are used in PSS method 
(m equal to 2). 

As shown in Figs. | and 2, the low order polynomial method 
(or uncorrected diffusion length method) is inadequate at short 
times for transient behavior. The method fails to converge to the 
exact solution at steady state. Further simulation shows that the 
solution of low order polynomial method will match the exact 
solution at long times only when the pore wall flux jn or 5(t) 
is a constant. Although the corrected diffusion length method 
provides improved accuracy at short times by incorporating the 
empirical exponential term, it has the same problem as the uncor- 
rected method at long times where the exponential term becomes 
close to zero (see Eqs. (6a) and (6b)). These simulation results 
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0.20 a High order polyn. 
Sa Low order polyn 
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Fig. 1. Comparison of approximation solutions with the exact solution of the 
dimensionless diffusion equation. 5(t) = 1 — T. 
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Fig. 2. Comparison of approximation solutions with the exact solution of the 
dimensionless diffusion equation. 5(t) = sin(5rT). 


indicate that the low order polynomial method and the corrected 
diffusion length method cannot provide as accurate results as the 
Duhamel’s superposition method, the PSS method or the high 
order polynomial method. 

By using a high order polynomial model for the concentration 
profile in the solid particles, the high order polynomial method 
is not only able to describe transient behavior in more detail, 
but also it is able to provide much more accurate results at long 
times. In addition, the method remains computational simplistic 
and can be easily implemented in a porous electrode model. 

The usage of the PSS method is mainly affected by the num- 
ber of summation terms included. If no summation terms are 
used, the method degrades to the low order polynomial method. 
Increasing the number of summation terms improves the accu- 
racy of the method, mostly at short times. Our simulation shows 
that PSS method with two or three summation terms is able 
to provide accurate results under a wide range of operating 
conditions. More summation terms require solving more differ- 
ential equations (Eq. (9b)) and could pose numerical difficulties 
because of the exponential terms. Fig. 3 shows that the PSS 
method with four summation terms must be solved with a tight- 
ened error tolerance in order to obtain accurate results when 
the dimensionless pore wall flux is 6(t)=sin(5rt). If the PSS 
method is used in the porous electrode model, the approxima- 
tion equations are tightly coupled to other equations for the liquid 
phase concentration and potentials. The smaller error tolerance 
required by the PSS method with several summation terms could 
have negative effects on solving the whole equation system, 
rendering it computationally more intensive. 

The Duhamel’s superposition method can provide the most 
accurate solutions when used with a small time step At. How- 
ever, this also leads to its disadvantages. The method requires 
significantly more computation time when the time frame is long 
and the time step is small. It also demands more storage space 
to record values of dependent variables from all previous time 
steps because those values are required in the calculation for the 
current time step (Eq. (2a)). When the method is used in porous 
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Cye-Cayg 


—— 4 Terms, Same error tolerance as 2 terms 
—<— 4 Terms, Tightened error tolerance 
—o— 2 Terms 


0.0 0.2 0.4 0.6 
T 


Fig. 3. A tight error tolerance is required to solve PSS method with more sum- 
mation terms. PSS method with four summation terms requires smaller error 
tolerance to get accurate results. ô(T) = sin(5T). 


electrode model and localized at each discretization node along 
spatial coordinates, the computation and storage requirement 
become exacerbated, as will be shown in the following section. 

The preliminary dimensionless analysis indicates that pseudo 
steady state method and the high order polynomial method 
would be good choices for the use in porous electrode model. 
In the following, we will compare the usage of the approxima- 
tion methods in a porous electrode model to predict discharge 
profiles of a LiCoO2 half cell. 


4. Approximation methods in a porous electrode model 


Equations for a porous electrode model have been described 
extensively in many papers [1,2] and they will not be repeated 
here. The discharge profiles of a LiCoO2 electrode were mea- 
sured in an experiment, the details of which were revealed 
elsewhere [18]. A pseudo-2D porous electrode model was used 
in the study to simulate the discharge profiles. By pseudo-2D, we 
mean that the model is solved in the spatial coordinate x along 
with the solid phase diffusion equation being solved in the extra 
r dimension at each x. In order to obtain good simulation results, 
the symmetry factor £ in the Butler Volmer equation was empir- 
ically assumed to change with the cell state of discharge [18]. 
The empirical expression for the Butler Volmer equation was 
formulated as follows: 


. OF Qa lc Oa F 
Jn = Kpce* (Cmax — Cse) “Coe (ex (Sr. — þe — Ueq 
in FR lad in FR 
— Jn £)) exp ( RT (bs — ge Ueq Jn D)) 
(12a) 
=ß=0.5 (: : ) (12b) 
comin 1 + exp(a(b — cse/Cmax)) 
a =1-8 (12c) 


Potential (V) 


0.4 0.6 0.8 1.0 
x in Li,CoO, 


Fig. 4. Experimental (symbols) and simulated (lines) discharge profiles at differ- 
ent current rates. The rates from top to bottom: C/7, C/2.7, and C/1.3 (C=4 mA). 
The pseudo-2D porous electrode model well captures the effect of increased 
solid phase diffusion limitation at higher current rates indicated by the increased 
curvature at the beginning of discharge. 


Fig. 4 presents a comparison of the experimental and 
simulated discharge profiles using the pseudo-2D porous elec- 
trode model at different current rates: C/7, C/2.7, and C/1.3 
(C=4mA). The pseudo-2D model captures well the effect of 
increased solid phase diffusion limitation at higher current rates, 
as indicated by the increased curvature at the beginning of dis- 
charge profiles. The parameters used in the simulation are listed 
in Table 1. 

The approximation methods are used next in the porous elec- 
trode model to eliminate the need to solve the solid phase 
diffusion equation in the r direction. These models use the 
same parameters given in Table 1 and the simulation results 
are compared with those obtained using the pseudo-2D model 
in Figs. 5-7. 


—o— Low order polyn. 
—o— Corrected diff. Lgth 
—+— High order polyn 
4. —— PSS 
30 —7— Duhamel 
—*— Pseudo2D 


— 4.25 
= 
S&S 
5 
5 4.20 
a 

4.15 

4.10 

0.35 0.40 0.45 0.50 0.55 
x in Li,CoO, 


Fig. 5. Comparison of simulation results from porous electrode models with 
approximation methods and pseudo-2D model. Discharge rate is C/7 (C=4 mA). 
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Table 1 
Model parameters used to simulate the discharge curves 
Parameter Value 
LiCoOy> electrode 

T (°C) 15? 

Wg) 0.0245? 

dp (cm) 64 x 10-44 

Ee 0.30? 

S (cm?) 1.267? 

R, (cm) 10 x 1074? 

os (Scm™!) 0.1? 

Brug 1.5° 

D; (cm? s7!) 1.2 x 10710¢ 

Rg (Q cm?) 200° 

X0,p 0.393° 

kp 2.59 x 1076° 
Separator 

Ce0 1 x 107-34 

ôs (cm) 25 x 1074? 

log(De) (cm? s7!) —4.43 —54/(T — 5 x 103 ce — 229) — 0.22 x 10? ce4 

Ke (S cm7!) Ce(—10.5 + 0.074T — 6.96 x 1075T? + 668ce — 17.8ceT + 0.028ceT? + 4.94 x 105c2 — 886c2 T)? 

v=0 -0 (14+ FE) 0.601 — 7.59c25 + 3.1 x 104(2.53 — 0.0052T)c!54 

S (cm?) 1.267? 

Ee 0.46" 

Parameters used in Eqs. (12a)—(12c) 
ae be 

CIT 21 0.92 
C/2.7 13 0.88 
C/1.3 8 0.8 


à Manufacturer data or experiment data. 
> From Refs. [3,19,20]. 

© Values fit to experiment data. 

4 From Ref. [21]. 


The simulation results show that the choice of an approxima- 
tion method is more important at high current rates than it is at 
low current rates. The long time behavior is also less affected 
than the short time behavior. The short time behavior is very 


—°— Low order polyn. 
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` —— Duhamel 
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Potential (V) 


“0.35 0.40 0.45 0.50 0.55 
x in Li,CoO, 


Fig. 6. Comparison of simulation results from porous electrode models with app- 
roximation methods and pseudo-2D model. Discharge rate is C/2.7 (C=4 mA). 


important in studying the discharge profiles of the LiCoO? elec- 
trode because it reflects valuable information about how severe 
the solid diffusion limitation is in the LiCoO> electrode. Sim- 
ulation results obtained from the model using the low order 


—— Low order polyn. 
—— Corrected diff. Lgth. 
—+— High order polyn. 
—o— PSS 
4.2 F —— Duhamel 
—— Pseudo2D 
= 
S 41b 
c 
g 
[e] 
oa 
4.0F 
3.9 
0.35 0.40 0.45 0.50 0.55 


xin Li CoO, 


Fig. 7. Comparison of simulation results from porous electrode models with app- 
roximation methods and pseudo-2D model. Discharge rate is C/1.3 (C =4 mA). 
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Table 2 
Computation time (in seconds) required to solve the porous electrode models* 
C/T C/2.7 C/1.3 
Low order polynomial 1 0.9 0.9 
Corrected diffusion length 1.1 1.0 1.0 
High order polynomial 1.2 1.1 1.1 
Duhamel 62 89 87 
PSS 1.5 1.4 1.3 
Pseudo-2D 10.5 10 9.6 


à Average of five runs. 


polynomial method (or uncorrected diffusion length method) 
have significant errors at the beginning of discharge, especially 
when the discharge current is high. The model basically failed 
to show how the solid phase diffusion influences the discharge 
curves of the LiCoO) electrode. The corrected diffusion length 
method provides significantly improved results compared to 
the uncorrected method in describing transient behaviors, but it 
tends to slightly over predict the cell potential. The model using 
the Duhamel superposition method is accurate, but it requires 
much more computation time than the pseudo-2D model. All 
other models using an approximation method have significant 
speed advantages over the pseudo-2D model. Table 2 lists the 
computation times required to solve these models using DASSL 
[22] solver. The model using the high order polynomial method 
provides good results at most times, but it is inferior to the PSS 
method at short times when the current is high. The PSS method 
stands out in the comparison because the model using the PSS 
method is able to provide simulation results not only with accu- 
racy but also with a speed advantage. It should be noted that even 
though the study is based on a LiCoO3 electrode with character- 
istic particles size of 10 x 1074 cm, the models are expected to 
perform mostly the same as long as the lithium transport in the 
particles can be mathematically described through Fick’s diffu- 
sion law (Eqs. (1a)—(1c)), regardless of the electrode materials 
and particle sizes. 


5. Summary 


Approximation methods are often used in porous electrode 
models to make them more computationally efficient. In this 
work, the Duhamel’s superposition method, the corrected diffu- 
sion length method, the polynomial approximation method, and 


the pseudo steady state method are compared to each other in 
a dimensionless analysis study, and they are used in a porous 
electrode model to describe the discharge profiles of a LiCoO2 
electrode in a half cell. The low order polynomial method was 
found to be inadequate for transient behavior in both studies. 
Although the corrected diffusion length method proved reason- 
ably accurate at short times, this method did not converge to 
the exact solution at steady state. The Duhamel’s superposi- 
tion method requires extensive computation time and storage 
space. Simulation results show that porous electrode model 
using Duhamel’s method requires even more computation time 
than the pseudo-2D model when the DASSL solver is used. 
This study indicates that the pseudo steady state method or the 
high order polynomial method should be considered first if an 
approximation method is to be used in a porous electrode model. 
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